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Abstract. An industrial scheme, to simulate the two compressible phase flow in porous 
media, consists in a finite volume method together with a phase-by-phase upstream scheme. 
The implicit finite volume scheme satisfies industrial constraints of robustness. We show 
that the proposed scheme satisfy the maximum principle for the saturation, a discrete energy 
estimate on the pressures and a function of the saturation that denote capillary terms. These 
stabilities results allow us to derive the convergence of a subsequence to a weak solution of 
the continuous equations as the size of the discretization tends to zero. The proof is given 
for the complete system when the density of the each phase depends on the own pressure. 



1. Introduction 

A rigorous mathematical study of a petroleum engineering schemes takes an important 
place in oil recovery engineering for production of hydrocarbons from petroleum reservoirs. 
This important problem renews the mathematical interest in the equations describing the 
multi-phase flows through porous media. The derivation of the mathematical equations de- 
scribing this phenomenon may be found in [6], [10]. The differential equations describing 
the flow of two incompressible, immiscible fluids in porous media have been studied in the 
past decades. Existence of weak solutions to these equations has been shown under various 
assumptions on physical data Ill[l0l[IIl[l2l[l3l[I71[l8l[2ll|^. 

The numerical discretization of the two-phase incompressible immiscible flows has been 
the object of several studies, the description of the numerical treatment by finite difference 
scheme may be found in the books [5], |27j . 

The finite volume methods have been proved to be well adapted to discretize conservative 
equations and have been used in industry because they are cheap, simple to code and robust. 
The porous media problems are one of the privileged field of applications. This success induced 
us to study and prove the mathematical convergence of a classical finite volume method for 
a model of two-phase flow in porous media. 

For the two-phase incompressible immiscible flows, the convergence of a cell-centered flnite 
volume scheme to a weak solution is studied in [26], and for a cell-centered flnite volume 
scheme, using a "phase by phase" upstream choice for computations of the fluxes have been 
treated in [16] and in The authors give an iterative method to calculate explicitly the 
phase by phase upwind scheme in the case where the flow is driven by gravitational forces 
and the capillary pressure is neglected. An introduction of the cell-centered flnite volume can 
be found in [15j. 

For the convergence analysis of an approximation to miscible fluid flows in porous media 
by combining mixed finite element and finite volume methods, we refer to [2], [3]. 
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Pioneers works have been done recently by C. Galusinski and M. Saad in a serie of articles 
about "Degenerate parabolic system for compressible, immiscible, two-phase flows in porous 
media" ([19], [20j, ^21j) when the densities depend on the global pressure , and by Z. Khalil 
and M. Saad in ([22], [23]) for the general case where the density of each phase depends on 
its own pressure. And for the two compressible, partially miscible flow in porous media, we 
refer to |9], [28]. For the convergence analysis of a finite volume scheme for a degenerate 
compressible and immiscible flow in porous media with the feature of global pressure, we 
refer to [7]. 

In this paper, we consider a two-phase flow model where the fluids are immiscible. The 
medium is saturated by a two compressible phase flows. The model is treated without sim- 
plified assumptions on the density of each phase, we consider that the density of each phase 
depends on its corresponding pressure. It is well known that equations arising from multi- 
phase flow in porous media are degenerated. The first type of degeneracy derives from the 
behavior of relative permeability of each phase which vanishes when his saturation goes to 
zero. The second type of degeneracy is due to the time derivative term when the saturation 
of each phase vanishes. 

This paper deals with construction and convergence analysis of a finite volume scheme for 
two compressible and immiscible flow in porous media without simplified assumptions on the 
state law of the density of each phase. 

The goal of this paper is to show that the approximate solution obtained with the proposed 



upwind finite volume scheme (3.8)-(3.9) converges as the mesh size tends to zero, to a solution 
of system (2.1) in an appropriate sense defined in section [2| In section [3| we introduce some 
notations for the finite volume method and we present our numerical scheme and the main 
theorem of convergence. 

In section |4j we derive three preliminary fundamental lemmas. In fact, we will see that we 
can't control the discrete gradient of pressure since the mobility of each phase vanishes in the 
region where the phase is missing. So we are going to use the feature of global pressure. We 
show that the control of velocities ensures the control of the global pressure and a dissipative 
term on saturation in the whole domain regardless of the presence or the disappearance of 
the phases. 

Section [5] is devoted to a maximum principle on saturation and a well posedness of the 
scheme which inspired from H.W. Alt, S. Luckhaus [Ij. Section [7] is devoted to a space-time 

compactness of sequences of approximate solutions. 
Finally, the passage to the limit on the scheme and convergence analysis are performed in 
section [H Some numerical results are stated in the last section [H 

2. Mathematical formulation of the continuous problem 

Let us state the physical model describing the immiscible displacement of two compressible 
fluids in porous media. Let T > be the flnal time flxed, and let be a bounded open subset 
of {i> 1). We set Qt = (0,T) x 0, = (0,T) x dn. The mass conservation of each 
phase is given in Qt 

(2.1) (l){x)dt{pa{Pa)Sa){t,x) + div{pa{Pa)^a){t, x) + Pa{Pa) Saf p{t, x) = Pa{pa)sifi{t,x), 

where 4>, pa and Sa are respectively the porosity of the medium, the density of the a phase 
and the saturation of the a phase. Here the functions fj and fp are respectively the injection 



and production terms. Note that in equation (2.1) the injection term is multiplied by a 



known saturation corresponding to the known injected fluid, whereas the production term 
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is multiplied by the unknown saturation Sa corresponding to the produced fluid. 
The velocity of each fluid Vq, is given by the Darcy law: 

(2.2) v„ = -K^^^^^(Vp„-p„(p,)g), a = l,g. 

where K is the permeability tensor of the porous medium, kr^ the relative permeability of the 
a phase, the constant a-phase's viscosity, pa the a-phase's pressure and g is the gravity 
term. Assuming that the phases occupy the whole pore space, the phase saturations satisfy 

(2.3) Sl+Sg = l. 

The curvature of the contact surface between the two fluids links the jump of pressure of 
the two phases to the saturation by the capillary pressure law in order to close the system 



(2.1)-(2.3) 



(2.4) Pc{si{t, x)) = pg{t, x) - pi{t, x). 

With the arbitrary choice of ( |2.4[ ) (the jump of pressure is a function of si), the application 
•5/ ^ Pc{si) is non-increasing, (^(s;) < 0, for all s; E [0,1]), and usually Pc{si = 1) = 
when the wetting fluid is at its maximum saturation. 

2.1. Assumptions and main result. The model is treated without simplified assumptions 
on the density of each phase, we consider that the density of each phase depends on its 
corresponding pressure. The main point is to handle a priori estimates on the approximate 
solution. The studied system represents two kinds of degeneracy: the degeneracy for evolution 
terms dt{paSa) and the degeneracy for dissipative terms div{paMaVpa) when the saturation 
vanishes. We will see in the section [5] that we can't control the discrete gradient of pressure 
since the mobility of each phase vanishes in the region where the phase is missing. So, we 
are going to use the feature of global pressure to obtain uniform estimates on the discrete 
gradient of the global pressure and the discrete gradient of the capillary term B (defined on 



(2.7)) to treat the degeneracy of this system. 

Let us summarize some useful notations in the sequel. We recall the conception of the 
global pressure as describe in [10] 



M{Sl)Vp = Ml{Sl)Vpi + Mg{Sg)Vpg, 

with the a-phase's mobility and the total mobility are defined by 

Ma{Sa) = krASa)/Pa, M{si) = Mi{si) + Mg{Sg). 

This global pressure p can be written as 

(2-5) P = Pg+p{si) =pi+p{si), 

or the artificial pressures are denoted by p and p defined by: 

(2.6) ^>') = -i, i4y^^^^^''""'^^''^=/o i4y"^^^^'^ 

We also define the capillary terms by 

.(.) = -«^f^(.)>o, 

M{si) dsi 
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and let us finally define the function B from [0, 1] to M by: 

W ^ r ( r Mi{z)Mg{z) dp,, 



Jo 



(2.7) = - / Mi(z)^^{z)dz = £ M,(z)^(z)dz. 







We complete the description of the model (|2.1l) by introducing boundary conditions and initial 



conditions. To the system (2.1)-(2.4) we add the following mixed boundary conditions. We 
consider the boundary dQ = F/ U Timp, where Ti denotes the water injection boundary and 
Timp the impervious one. 

^^g^ r pi{t, x)=pg{t,x) = on {0,T)xTi, 

I piVl ■ n = pgYg • n = on (0, T) X Timp, 

where n is the outward normal to Timp- 

The initial conditions are defined on pressures 

(2.9) pa{t = 0)=Pa for a = l,g inn 

We are going to construct a finite volume scheme on orthogonal admissible mesh, we treat 
here the case where 

K = kla 

where A; is a constant positive. For clarity, we take k = 1 which equivalent to change the scale 
in time. 

Next we introduce some physically relevant assumptions on the coefficients of the system. 

(HI) There is two positive constants 0o ^-nd (pi such that < (j){x) < (pi almost everywhere 
X £n. 

{H2) The functions Mi and Mg belongs to C°([0, 1],M+), Ma(s„ = 0) = 0. In addition, there 
is a positive constant mo > such that for all si G [0, 1], 

Mi{si) + Mg{sg) > mo. 

{H3) {fp,fj) G {L'^{QT)f, fp{t,x), fj{t,x) > almost everywhere {t,x) £ Qt- 

(i74) The density pa is C^(M), increasing and there exist two positive constants pm > and 

PM > such that < Pm < Pa{Pg) < PM- 

{H5) The capillary pressure fonction pc G C^([0, 1];M+), decreasing and there exists pc > 

such that <pc< \ ^\- 
{H6) The function 7 G ([0, 1]; R+) satisfies 7(5/) > for < < 1 and 7(5^ = 1) = 

7(5; = 0) = 0. We assume that B~^ (the inverse of B{si) = J^' 7(2:)d2:) is an Holdeij^ 

function of order 6, with < 9 < 1, on [0,^6(1)]. 

The assumptions (lQ-(]^6]) are classical for porous media. Note that, due to the boundedness 



of the capillary pressure function, the functions p and p defined in (2.6) are bounded on [0, 1]. 
Let us define the following Sobolev space 

H^^{n) = {n G H\n);u = sur Ti}, 

this is an Hilbert space with the norm ||n||j|^i (j^) = ||Vn||(2,2(f7))« . 



-'^This means that there exists a positive constant c such that for all a, 6 G [0, one has \B ^(c 

B-\b)\ < c\a-bf. 
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Definition 1. [Weak solutions) . Under assumptions and definitions (2.5)-(2.9) 

with the fact that p^, p^ belongs to L'^{Q) and satisfies < s^ < 1 almost everywhere in 
Q, then the pair {pi,Pg) is a weak solution of problem ( |2.ip satisfying : 

(2.10) Pa G L\0, T; L^{n)), y^M^Vpa € (L^O, T; L\n))Y, 

(2.11) < sait,x) < 1 a.e in Qt (a = l,g), B{si) G L'^{0,T; H^^{n)), 

(2.12) ^dt{pa{Pa)sa) G L\0,T; (H^^in))') e L^{0,T;{H^^{n)y), 
such that for all ip, e C^{[0,T]; H^^{n)) withip{T,-) = ^p{T, •) = 0, 

(t)pi{pi)sidt(pdxdt - / 4>{x)pi{p^{x))s^{x)(p{0,x)dx 

Qt 

(2.13) + / Mi{si)pi{pi)Vpi-V^dxdt- I Mi{si)pf{pi)g-V^dxdt 

JQt JQt 

+ / pi{pi)sifp(pdxdt= / pi{pi)slfj(fdxdt, 
JQt JQt 

- I (t)Pg{Pg)Sgdttpdxdt - / ^ (x ) (^^ (x ) ) (x) ( , X ) rfx 

JQt J^ 

(2.14) + / Mg{Sg)Pg{Pg)Vpg -Vlljdxdt - [ M g {s g) P^Pg) g ' V ^^dxdt 

JQt JQt 

3. The finite volume scheme 

3.1. Finite volume definitions and notations. Following |15j . let us define a finite volume 
discretization of x (0, T). 

Definition 2. {Admissible mesh ofO,). An admissible mesh TofQ is given by a set of open 
bounded polygonal convex subsets of 17 called control volumes and a family of points (the 
"centers" of control volumes) satisfying the following properties: 

(1) The closure of the union of all control volumes is fi. We denote by \K\ the measure 
of K , and define 

h = size{T) = max{diam{K) , K £ T}. 

(2) For any {K, L) G with K ^ L, then Kr\L = %. One denotes by £ CT^ the set 
of {K, L) such that the d — 1-Lebesgue measure of K D L is positive. For (K, L) € 8, 
one denotes (7k\l = KCiL and the d— 1-Lebesgue measure of ax\L- ^"^d one 
denotes rj^^i the unit normal vector to ctj^^^ outward to K 

(3) For any K €z T, one defines N{K) = {L £ T,{K,L) G £} and one assumes that 
dK = K\K = (KndQ) U {ULeN(K)^K\L)- 

(4) The family of points {xk)k<^t is such that xk G K {for all K G T) and, if L £ 
N{K), it is assumed that the straight line {xk,xl) is orthogonal to ctx\l- We set 

dx\L = d{xK,XL) the distance between the points xk and xl, and Tx\l = ^z^^? ^^'^^ 
is sometimes called the "transmissivity" through (Jx\l (see Figure^. 
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(5) Let ^ > 0. We assume the following regularity of the mesh : 

(3.1) Vi^er, WK\L\dK\L<m 

L&N{K) 




Figure 1 . Control volumes, centers and diamonds 

We denote by Hh{Q) C L'^{Q) the space of functions which are piecewise constant on each 
control volume K £ T. For all Uh € Hh{^) and for all K £ T, we denote by uk the constant 
value of Ufi in K. For {uh,Vfi) S {Hh{Q))'^, we define the following inner product: 

{uh,Vh)H^ = ^ "Y Y ^—7^{UL - Uk){vl - Vk), 
K€TLeN{K) 

and the norm in Hh{Q) by 

Finally, we define Lfi{^) C L'^{Q,) the space of functions which are piecewise constant on each 
control volume K £ T with the associated norm 

iuh,Vh)L,^(^n) = \K\UKVK, = Y \K\\uK? , 

Ker KeT 

for {ufi,Vh) £ Further, a diamond T^|i is constructed upon the interface fT^i^, 

having xk, xl for vertices (see Figure [l| and the ^-dimensional mesure of Ti^\i equals 

to J \(yK\L\ dK\L- 
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The discrete gradient V/^n/j of a constant per control volume function is defined as the 
constant per diamond T^i^, M^-valued function with values 



And the semi-norm coincides with the L'^{Q) norm of V ^Uh, in fact 

KeT LeN{K)-^'^K\i^ KeTLeN{K) ^ 



'''K\L\ 



ESr^ ]^K\l\. |2 II ||2 



L&N(K) 



d 



K\L 



We assimilate a discrete field o^i ^ to the piecewise constant vector-function 

The discrete divergence of the field is defined as the discrete function Wh = divhFh with 
the entires 

(3.2) divKFh:=i^^ Y \^k\l\^k\l ■ Vk\l- 

LeN{K) 

The problem under consideration is time-dependent, hence we also need to discretize the 
time interval (0, T). 

Definition 3. (Time discretization). A time discretization o/(0,T) is given by an integer 
value N and by a strictly increasing sequence of real values (i")ne[o,v+i] with t^ = and 
t^^^ = T. Without restriction, we consider a uniform step time 6t = t"''^^ — i", for n € [0, N]. 

We may then define a discretization of the whole domain Q x (0, T) in the following way: 

Definition 4. (Discretization of Q x (0,T)). A finite volume discretization T> of 0, x (0,T) 
is defined by 

V = (t, £, (xK)Ker, N, (i")„e[o,v]) , 

where F,£, (xk)k£T '^i^ admissible mesh ofQ in the sense of Definition^and N, (i")ne[o,Af] 
is a time discretization of (0,T) in the sense of Definition^ One then sets 

size(T>) = max(size(T) , St) . 

Definition 5. (Discrete functions and notations). LetT> be a discretization of x (0,T) in 
the sense of Definition^ We denote any function from T x [0, A'^ -|- 1] to M by using the sub- 
scripts, (sa,i) and Pa,!) for instance) and we denote its value at the point (xK,t^) using the 
subscript K and the superscript n (s^ ^ for instance, we then denote Sa,v = (^a K)KeT,ne[o,N+i]] 
To any discrete function ux> corresponds an approximate function defined almost everywhere 
onQx (0,T) by: 

uv(t,x) =n^+\ for a.e. (t,x) G (t",t"+^) x K,\/K e T,'in G [0,A^]. 
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For any continuous function / : M — t- M, /{uj}) denotes the discrete function {K,n 
f{u^^). if L G N{K), and uj) is a discrete function, we denote by = — 

For example, 5l+l{f{u)) = f{ul+') - f{ul+'^ 



Let us recall the following two lemmas : 

Lemma 1. (Discrete Poincare inequality ) [15j. Let Q be an open bounded polygonal subset 
of^^, 1 = 2 or 3. Let T be a finite volume discretization of in the sense of Definition^ 
and let u be a function which is constant on each cell K £ T, that is, u{x) = uk if x £ K, 
then 

\\u\\L^(n) ^ diam{n) \\u\\h^^^) , 
where |M|//^(f7) is the discrete Hq norm. 

Remark 1. {Dirichlet condition on part of the boundary). The lemma ^ gives a discrete 
Poincare inequality for Dirichlet boundary conditions on the boundary d^l. In the case of 
Dirichlet condition on part of the boundary only, it is still possible to prove a discrete Poincare 
inequality provided that the polygonal bounded open set $7 is connected. 

Lemma 2. (Discrete integration by parts formula). Let Fj^/^, K £ T and L € N{K) be a 
value in M depends on K and L such that Fj^/j^ = —Fj^^j^ and let ip be a function which is 
constant on each cell K £ T, that is, ip{x) = ipx if x £ K, then 

(3-3) Fk/l^k = -]^Y 5Z ^K/iifL-fK) 

Ker LeN{K) KeT LeN{K) 

Consequently, if Fj^jj^ = ax/ii^L — bx), with a^ji^ = aijx, then 

X] X] "'K/L{bL-bK)'^K = -]^Y ^ aK/LibL-bK){^L-fK) 
KeT L&N{K) KeT LeN{K) 

3.2. The coupled finite volume scheme. The finite volume scheme is obtained by writing 
the balance equations of the fluxes on each control volume. Let P be a discretization of 
Q X (0, r) in the sense of Definition |4j Let us integrate equations (2.1) over each control 



volume K. By using the Green formula, if $ is a vector field, the integral of div(<^) on a 



control volume K is equal to the sum of the normal fluxes of $ on the edges (3.2). Here 
we apply this formula to approximate Ma{sa)^Pa ■ flK\L-, (o^ = ^5) by means of the values 
Sa,K,Sa,L and Pa,K,Pa,L that are available in the neighborhood of the interface cr^^i. To do 
this, let us use some function Ga of (a, b, c) G M'^ . The numerical convection flux functions 
Ga £ C(M^,M), are required to satisfy the properties: 

(a) b, c) is non-decreasing for all b,c £M, 

and Ga{a, •, c) is non-increasing for all a, c G M; 
(3.5) { (b) Gaia, a, c) = -Ma{a) c for all a, c G M; 

(c) Ga{a, b, c) = —Ga{b, a, — c) and there exists C > such that 
\Gaia,b,c)\ < C {\a\ + \b\)\c\ for ah a, 6, c G M. 

Note that the assumptions (a), (b) and (c) are standard and they respectively ensure the 
maximum principle on saturation, the consistency of the numerical flux and the conservation 
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of the numerical flux on each interface. Practical examples of numerical convective flux 
functions can be found in [T5l. 



In our context, we consider an upwind scheme, the numerical flux Ga satisfying (3.5) 
defined by 

(3.6) Gc,(a,6,c) = -Mc,(6)c+ + Mc,(a)c" 

where c"*" = max(c, 0) and c~ = max(— c, 0). Note that the function Sa ^ Ma{sa) is non- 
decreasing, which lead to the monotony property of the function Ga- 

The resulting equation is discretized with a implicit Euler scheme in time; the normal 
gradients are discretized with a centered finite difference scheme. 

Denote by Pa,v = {p^K)KeT,ne[G,N] and s„,x, = (s^^)i^er,ne[o,Ar] the discrete unknowns 
corresponding to pa and Sq,. The finite volume scheme is the following set of equations : 



(3.7) 



pI,k = ^ j^pl{x)dx, 8% = ^\j^ sl{x)dx, for aU K G T, 



(3.8) \K\(pK + 2^ T^K\LPi^K\L^ly^l,K '^l,L '^K\lW)) 

L(^N{K) 

+ Kk' + 1^1 piip?;KKi'f^s = \K\ Pi{p?;MKr^'f?;K\ 

(3.yj \K\(I)K h 2^ T^K\LPg^K\L^9\^g,K^^g,L ^°K\LyP9)) 

LeN(K) 

+ F,:^' + 1^1 pMS^TKfp^K = \K\ PApTK)i4,Kr^'f?j<'^ 

(3.10) P.(^^S)=P"i'-P,"S, 

where F-^ (a = I, ,) the approximation of / p^(K-)M„(.r^)g • ,,„dr(x) by an 

JdK 

upwind scheme: 



(3.11) i^:^= E OV= E i^i^iLi(CiiL)%^-(CK)gxiL-^«(CL)gi.iL 

LeN(K) LeN{K) 

with g^i^ := (g • Vkil)^ and g^j^;^ := (g • riK\L) - Notice that the source terms are, for 
n e {0, . . . , iV - 1} 



The mean value of the density of each phase on interfaces is not classical since it is given 



as 



(3.12) 



n+l 
Pa,K\L 



1 rPa,L 

n + l_„n + l J„n + 1 



Pc,L 
1 

n+l 



Pc{C) 



dc -^^pi:k^p 

otherwise. 



n+l 
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This choice is crucial to obtain estimates on discrete pressures. 

Note that the numerical fluxes to approach the gravity terms Fq, are nondecreasing with 
respect to Sa.K and nonincreasing with respect to Sa,L- 

The upwind fluxes (3.6) can be rewritten in the equivalent form 

(3.13) G^is-^lsi;-,';5lY,{p^)) = -M„(.^5| J S^^Y^iPa), 

where Mq,(s^^|^) denote the upwind discretization of Ma{sa) on the interface a^^j^ and 

^ ^ otherwise, 
with the set £a~^^ subset of S such that 

(3.15) £^+' = m L) G £, = pl-^l - vl^^ < 0}. 

We extend the mobility functions Sq, i— )• Ma{sa) outside [0, 1] by continuous constant functions. 
We show below (see Prop. [2]) that there exists at least one solution to this scheme. From this 
discrete solution, we build an approximation solution pa^v defined almost everywhere on Qt 
by (see Definition [s]) : 

(3.16) PaMt,x)=PatK, e e (r,r+^). 

The main result of this paper is the following theorem. 

Theorem 1. Assume hypothesis hold. Let {2?m}mGN o, sequence of discretiza- 

tion of Qt in the sense of definition^ such that livum^+oo size{Vm) = 0. Let (Pa,s[^) G 
L^(Q,M) X L°°(r2,M). Then there exists an approximate solutions {pa,T>^)mGN corresponding 
to the system ^3^-{3.9), which converges (up to a subsequence) to a weak solution pa of 
(2.1) in the sense of the DefinitionUi 



4. Preliminary fundamental lemmas 

The mobility of each phase vanishes in the region where the phase is missing. Therefore, 
if we control the quantities MqVpq, in the L^-norm, this does not permit the control of the 
gradient of pressure of each phase. In the continuous case, we have the following relationship 
between the global pressure, capillary pressure and the pressure of each phase 

(4.1) M|Vpp + ^^\Vpc\' = M;|Vpz|2 + Mg\Vpf. 

This relationship, means that, the control of the velocities ensures the control of the global 
pressure and the capillary terms B in the whole domain regardless of the presence or the disap- 
pearance of the phases. This estimates (of the global pressure and the capillary terms B) has 
a major role in the analysis, to treat the degeneracy of the dissipative terms dW{paMa'Vpa) ■ 

In the discrete case, these relationship, are not obtained in a straightforward way. This 
equality is replaced by three discrete inequalities which we state in the following three lemmas. 

We derive in the next lemma the preliminary step to proof the estimates of the global 
pressure and the capillary terms given in Proposition [T] and Corollary [T] These lemmas are 
first used to prove a compactness lemma and then used for the convergence result. 
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Lemma 3. (Total mobility and global pressure [16jj. Under the assumptions (-H[T]) — (-H|6]) 
and the notations (2.5). Let T> be a finite volume discretization o/fi x (0,r) in the sense of 
Definition^ Then for all {K, L) ^ £ and for all n G [0, N] the following inequalities hold: 



(4.2) M-^, + M;;^\, > mo, 

and 



(4-3) -o(5^2(p))' < KM^K]liPi)f + 

The proof of this lemma is made by R. Eymard and al. in [16j. The proof of this resuh can 
be apphed for compressible flow since the proof use only the definition of the global pressure. 

Lemma 4. (Capillary term B). Under the assumptions (^/[T]) — (f/[6]) and the notations (2.5). 
LetV be a finite volume discretization of Q, x (0, T) in the sense of Definition^ Then there 
exists a constant C > such that for all (K, L) G £ and n £ [0, A'"] .• 

(4.4) {5-^Y,iB{s,))f < MpJ^l^Sl^lip^y + M;;^^],{6-^1{p, 

In the incompressible case (see [16]) this kind of estimate is obtained by using the mass 
conservation equation and under hypotheses ont the relative permeability of the a phase, 
whereas, the compressibility add more difficulties, our approach use only the definition of 
the function B and consequently this lemma can be used for compressible and incompressible 
degenerate flows. 

Proof. We take the same decomposition of the interface as that proposed by R. Eymard and 
al. in [m, namely the different possible cases {K,L) £ £l'+^n£^+^, {K,L) ^ £p+^U£^+^, 
{K, L) £ £p+'^ and (K, L) i £^+\ and the last case (K, L) ^ £1'"^^ and (K, L) E £^+^; where 
the sets £J'^^ and £-^+^ are defined in (|3.15l). We establish for the four cases. 

•First case. If {K,L) ^ £i and {K,L) £ £g. We may notice that if the upwind choice is 
different for the two equations, we have 

M"+.V =, max M„. 



a,K\L 



By definition of B in (2.7), there exists some a £ [si^k, si^l] such that 
we then get 



< 



•Second case: The case (K, L) £ £i and [K, L) ^ £g is similar. 
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•Third case: The case {K,L) £ £i and {K,L) £ Eg. We have 

(4.5) = (m,(.I:+^) + M,(.«J)) ))' 

- 2M,(.^J)5^2(^^)^^|2(P(^/)) - 2M,(si:+i)<5^2(p)<5^1^^(p(sO). 
We win distinguish the case s]'~^^ < s"^^ and the case > s"^^ 

(1) If we assume that < s"^^, we deduce that 

(a) S^l{p{si)) < since p{si) is nonincreasing, 

(b) '5^j'2(p('5/)) > since is nondecreasing, 

(c) = + S^kTl^P^si)) < 0. 
One then gets from (14. 5|) that: 



(4.6) 



>^MK.i:j^)+M,(.;^^)j(5-^^(p))^ 
+ M,(.;+i)(5^2(^5(^o)') + M,(.i:+^)(<5^i^^(p(.o))' 

-2MKsI:+^)5^2(p)'^k|l'(^^(«0)- 
The previous inequahty gives: 

+ M,(.^J)(5^2(^5(«/)))' + MK.[:i^)(5^2(p(^0))' 

+ Mi{sl+'W^+l{p)f + Mi^^^^^^^ 
which implies the inequahty: 

Or, by definition of B ( |2.7[ ), there exists some a G si^l] such that (5^+2 ('^(««)) 
^9(«)'^^t2(^(^0), we get then 



n+l 



< Mi{s?^'){6-]l{p0? + M,«+i)(5^+nP.))'- 



which is Q in that case. 
(2) If we assume that s"^^ < s"^^, we get that 



(4.7) 
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(a) S^^{p{si)) > since p{si) is nonincreasing, 

(b) S^^{p{si)) < since p{si) is nondecreasing, 

(c) dlYjp) = ^ViliPa) + ^K^liPisi)) < 0. 
One then gets from ( |4.5[ ) that: 

+ M,(.^J)(5^2(^5(^0)') + MK<+^)(5^2(^^(^0))' 
The previous inequahty gives: 

+ M,(.^+i)(<5^2(^5(^0)') + M,(sj;+i)(<5^1^^(p-(s0))2 

< ms?;K'){6Vilipi)? + M,(.^y)(5^2(p.))' 

+ Mgis-g^K^K^LiP)? + M,(.;+^)(5^1^^^5(./)))^ 

which imphes the inequahty: 

< + M,(.^S)(<5^2(^'.))'- 

Or, by definition of B (2.7) there exists some a S [si^i^jS/^^] such that 6^^{B{si)) = 
-Mi{a)6'p^l{p{si)), we get then 

< Mi{sli'){6-^l{p0)' + A^,(.^S)(5^2(p,))^ 

which is (j4]) in that case. 

•Fourth case: The case (K, L) ^ £i and {K, L) ^ £g is similar of the third case. □ 

Lemma 5. {Dissipative terms). Under the assumptions — and the notations (2.5). 
LetV be a finite volume discretization of fix (0, T) in the sense of Definition^ Then there 
exists a constant C > such that for all (K, L) £ £ and n € [0, N] 

(4-8) M-^l{S-^^l{p{s,))f < M5f45^2(P^))' + ^^SlLfeH^^))'' 

and 

(4-9) M;;^^],{6-^Y,ip{s,))f < M-^l[5],f,{p,)y + M;;l,],[6-YMy ■ 

Proof. In order to prove (4.8 ) and (4.9 ), we consider the exclusive cases {K, L) € <S"^^ n<S^+^, 
{K, L) ^ U (K, L) i 8'^+^ and {K, L) € f^+i and the last case {K, L) G £^+^ and 
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First case. If {K, L) ^ £i and {K, L) & Eg. We have 

M"t'V = max Ma, 

' [si.k,si,l\ 

and by definition of^there exists some a £ [sik, si,l] such that S'^lipisi)) = Mg {^)+Mi (a) ^k\l (^0 ) ; 
we get then 



which gives (4.8). For the discrete estimate (4.9) and by definition of p there exists some 
b G [si,k,si^l] such that d'^+liPisi)) = - (i^+S^) S'^l {Pc{si)), we get then 



which gives (4.9). 

Second case. The case {K,L) € £i and {K,L) ^ E'g is similar. 

The third case and the fourth case can be treated as the cases in the lemma [H □ 

5. A PRIORI ESTIMATES AND EXISTENCE OF THE APPROXIMATE SOLUTION 

We derive new energy estimates on the discrete velocities Ma(s^^|^)(5^j'2(pa). Never- 
theless, these estimates are degenerate in the sense that they do not permit the control of 
S^liPa), especially when a phase is missing. So, the global pressure has a major role in the 

analysis, we will show that the control of the discrete velocities Ma{s^'^j^^^)6^^{pa) ensures 
the control of the discrete gradient of the global pressure and the discrete gradient of the 
capillary term B in the whole domain regardless of the presence or the disappearance of the 
phases. 

The following section gives us some necessary energy estimates to prove the theorem [T] 

5.1. The maximum principle. Let us show in the following Lemma that the phase by 
phase upstream choice yields the stability of the scheme which is a basis to the analysis 
that we are going to perform. 

Lemma 6. {Maximum principe). Under assumptions Let {s^ j^)KeT ^ [0; 1] 

and let V = [T,£,{xK)KeT^^A't^)nelo,N]) ^e a discretization of Q x (0,T) in the sense of 



Definition^ and assume that {pa,v) is a solution of the finite volume (3.7)-(3.10). Then, the 
saturation [s^ j^)xeT,ne{o,...,N} remains in [0,1]. 

Proof. Let us show by induction in n that for all S T, > where a = l,g. For 

a = I, the claim is true for n = and for all K € T- We argue by induction that for 
all K G T, the claim is true up to order n. We consider the control volume K such that 
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s^^+i = min{sj^+^} and we seek that s[^+^ > 



For the above mentioned purpose, multiply the equation in (3.8) by — (sfj^ ) , we obtain 



(.1, 

L&N(K) 

- \K\ piiP^KKX^fpS^^K)- = - \K\ PiiP?;K')i4,Kr^'f?:'KiCKr < o- 

The numerical flux Gi is nonincreasing with respect to s"j^^ (see (a) in ( |3.5| )), and consistence 



(see (c) in (3.5)), we get 



(5-2) = -S^^YiiPi) Miisl^') (sl^'r = 0. 

Using the identity s"J^ = {s]^x^)^ — (sI^k^)^ , and the mobility M/ extended by zero on 
] - oo, 0], then Mz(s"+^)(s"^^)- = and ' 

(5.3) - f(^+^)(.^+1)- - \K\ piiPtK^X^rpSi^K)- 

L&N{K) 



Then, we deduce from (5.1) that 

and from the nonnegativity of s^xi obtain (s"x^)~ = 0. This implies that > and 

< < slj;^ for all n£[0,N-l] and LgT. 
In the same way, we prove Sg~^ > 0. □ 
5.2. Estimations on the pressures. 



Proposition 1. Let pa^v be a solution of (3.7)-(3.10). Then, there exists a constant C > 0, 
which only depends on Ma, ^, T , p^^, s^, s^, fp, fi and not on V, such that the following 
discrete L'^{0,T] H^{Q,)) estimates hold: 

N-l 

(5-4) E'^^E E ^K\LMa{s-^^^,)K-:2 -plSi' <c, 

n=0 KeT LeN{K) 

and 

N-l 

(5.5) E^^E E -k\l\pI-''-pT'?<C. 

n=0 KeT LeN{K) 
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Proof. We define the function TiaiPa) '■= Pa{Pa)ga{jPa) — Pa-, 'Pc{si) '■= f^^ Pc{z)dz and 

Pa (2) 



ga{Pa) = /(f" iz') '^-^- -'^'^ following proof, we denote by Cj various real values which 



only depend on Ma, ft, T, p^, s^, s^, fp, fi and not on V. To prove the estimate (5.4), we 
multiply (3.8) and (3.9) respectively by gi{pi,K), gg{Pg,K) and adding them, then summing 
the resulting equation over K and n. We thus get: 

(5.6) Ei + E2 + E^ + Ei = 0, 
where 

N-l 

^^=Y.Y. 1^1 ^K{{pi{ptK)stK - PliP?,K)s?,K) gi{ptK) 
n=0 K&T 

N-l 

n=0 KeT LeN{K) 
N-l 

n=0 KeT 
N-l 

n=0 A'gT 

To handle the first term of the equality (5.6). Let us forget the exponent n + 1 and let note 
with the exponent * the physical quantities at time t". In [22] the authors prove that : for 
all Sa >0 and s* > such that si + Sg = s* + s*g = 1, 

(5.7) {piipi)si - pl{pt)st)giipi) + {pg{Pg)Sg - Pgip*g)s*)ggiPg) 

> ni{pi)si - ni{pt)st + ng{pg)sg - ng{p;)s; - vm) + vM). 

The proof of (|5.7l) is based on the concavity property of Qa and Vc- So, this yields to 



(5.8) E,>Y,<^K \K\ [sflMpf^K) - sIk^pIk) + sIMpIk) " sIMpIk) 
KeT 

-Y.^K \K\ V.isf^j,) +Y.^K \K\ V.islj,). 
KeT KeT 

Using the fact that the numerical fluxes Gi and Gg are conservative in the sense of (c) in 
(3.5), we obtain by discrete integration by parts (see Lemma [2]) 

N-l 

^2 = 2 E E E ^K\L (<jtG'K.i:+\ ; - miPU')) 

n=0 KeT LeN{K) 
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and due to the correct choice of the density of the phase a on each interface, 

(5.9) ptmia.ivtk) - UpI^D) = p^K - <L , 

we obtain 

N-l 

^2 = 9E^*E E 



=0 KeTLeN{K) 



+ Gg{s 



n+l n+1. ^n+l 



The definition of the upwind fluxes in (3.13) imphes 

- pU') + G,{s-^^\s-X'-,6-Y,{pMpTk - pTl) 

Then, we obtain the fohowing equahty 

N-l 

(5-10) E, = -Y,StY: E + 

71=0 KeT LeN(K) 



To handle the other terms of the equality (5.6), firstly let us remark that the numerical fluxes 
of gravity term are conservative which satisfy F^^^j^ = —FJ^^]^ and -P^^^Il ~ ~'EgL^K^ ^° 
integrate by parts and we obtain 

Af-l 

^3 = ^ E E E WKiL\{F!;;;^l\9i{prK')-9^^^^^^^^ 

n=0 KgTL£N{K) 



According to the choice of the density of the phase a on each interface (5.9) and the definition 
(3.11 ) we obtain 

N-l 

^3 = -2E'^*E E kxiz.l/^SiL[^/Ki^)g;.|,-^/Ki^)gxiL]('^a^'(^'0) 

n=0 KeT LeN(K) 

1 ^-^ 

-2E'^*E E \-KiL\p';i]Lms-:ji)gj,\,-M,^^^^^^^^ 

n=0 KeT LeN{K) 

Recall the truncations of ^^liPa) 



with 6l+lipa) = (<5^2(Pa))+ - {Sl]l{Pa))-. So wc obtain 

7V-1 

n=0 KeT L&N{K) 
N-l 

+ 2E'^*E E ki^i^l/'SiL^'Kr)g..|L(^^2(pO)+ 

n=0 K&T LeN(K) 
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N-1 

2 



n=0 A'eT LeN{K) 
N-1 



2 

n=0 A'eT LeN{K) 

II 1 - 

From the following equality = {dx\LWK\L\)^TK\L ^^'^ ^-PP^Y the Cauchy-Schwarz in- 

equality to obtain 

N-1 

E3<2CY,StY, E dK\LWK\L\ 
n=0 ATer LeN{K) 

N-1 

n=0 XeT L€N(K) 

< 2CT\n\ 

n=0 K&T L&N{K) 

From the definition of the truncations of 5^^{pa), we obtain 

N-1 

(5.11) E,<2CT\n\ + -Y.^^Yl E 

n=0 ATer L(iN{K) 



+ ^.(CA-|L)('^^|L'(f.))')- 



The last term will be absorbed by the terms on pressures from the estimate (5.10). 

In order to estimate E^, using the fact that the densities are bounded and the map ga is 

sublinear {a..e.\g{pa)\ < C'IPal), we have 



N-1 



n=0 KeT 

then 

N-1 

\E,\ <CrY.6tY^ \K\ + /,';+i)(2|p^+i| + + 

ri=0 ATer 

Hence, by the Holder inequality, we get that 

N-1 

\E,\ < c, ii/p + ML.(Q^) ( 5t iIpM'.(^) ) ^ 

n=0 
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and, from the discrete Poincare inequality lemma [T| we get 

N-l 

(5.12) \E^\ < C73( ^ 5t IIpMh, y+CA- 

n=0 

The equality (5.6) with the inequalities ( |5.8| ), ( |5.10 ), (5.11), (5.12) give (5.4). Then we deduce 
( [5^ from □ 

We now state the following corollary, which is essential for the compactness and limit study. 

Corollary 1. From the previous Proposition, we deduce the following estimations: 



(5.13) 

(5.14) 

and 
(5.15) 



N-l 

E'^^E E rK\LKYL^13{si))f<C, 

n=0 KgT L&N{K) 
N-l 

E^*E E -/^il^Sil(^^2(^^(^o))^<^^> 

n=0 Xer L&N{K) 
N-l 

n=0 K&r LgN{K) 



□ 



Proof. The prove of the estimates (5.13), (5.14) and (5.15) are a direct consequence of the 
inequahty (Q, ([48|, (|49]) and the Proposition [Ij 

6. Existence of the finite volume scheme 



We start with a technical assertion to characterize the zeros of a vector field which stated 
and proved in 



Lemma 7. f|14j . p. 529) Assume the continuous function v : 

v{z) • z > if \\z\\ = r, 
for some r > 0. Then there exists a point z with \\z\\ < r such that 

v{z) = 0. 



satisfies 



Proposition 2. The problem (3.8)-(3.9) admits at least one solution {p]^x^Pg K){K,n)e'D- 
Proof. At the beginning of the proof, we set the following notations; 

M := Card{T), 

PlM := Kinder GM^, 



We define the map Th '■ 

%{piM'Pg,M) = {{'Ti,K}KeT,{'Tg,K}KeT) where 
pAp^kKk - P^^pIkKk 



Ti. 



LK 



\K\ 



6t 



LeN{K) 
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PgiP. 



n+l\ n+l 



Pg{P\ 



'g,K>^g,K 



6t 



I \ ^ _ n+l / n+l n+l xn+l / \\ 
+ 2^ ^K\LPg^K\LGg{Sg^K^S9,L^^K\LiPg)) 
LeN{K) 



+ FgT + 1^1 PsiPg:K)iCKfp:K - (4,X)"+V,"i^). 



Note that Th is well defined as a continuous function. Also we define the following homeo- 
morphism : 



iM X M.M ^^^^ ^Yiat. 



^{piMiPgM) = {viM^'^gM) 

where Va,M = {5a(K^)}ii-er- 

Now let us consider the following continuous mapping Vh defined as 

Vh{vi,M,Vg,M) = Tho J^'^ivi^MiVg^M) = Th{piM^PgM)- 
According to Lemma [7j our goal now is to show that 

(6.3) Vh{vi,M,Vg^M) ■ {viM^'^gM) > ^oi \\{vi^M,Vg^M)\\M?M = r > 

and for a sufficiently large r. 
We observe that 

Vh{viM^Vg,M) ■ {viM^VgM) >^ E '^^ 1^1 fe'^Ki') - <KnPY,K) 



sg^'np';^') - ^unPiK: 



6t 



n+i 11" ^ 



for some constants C > 0. This implies that 

Vhivi,M,Vg,M) ■ ivi,M,Vg,M) ^ " ^ E '^^ '^1 {<K^(pIk) + sIk^PIk)) 



(6.4) 



5t 



for some constants C, C > 0. Finally using the fact that is a Lipschitz function, then 
there exists a constant C > such that 



a9i{p?,i')}Ker,{9g{pTK'))}Ker) 



< C 



PLh 



+ 



Pg,h 



L2(Q) 



-'^'-^yWPh \\L2(n)^ \\Ph \\L'2{n)^ \\Ph \\L'2{n) 



Using this to deduce from (6.4) that (6.3) holds for r large enough. Hence, we obtain the 
existence of at least one solution to the scheme (3.8)-(3.9). □ 
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7. Compactness properties 

In this section we derive estimates on differences of space and time translates of the func- 
tion (j)vPa{Pav)saV which imply that the sequence (j)vPa{Pav)saV is relatively compact in 

lHQt). 

We replace the study of discrete functions Ua,v = 4>vPa{'Pa,v)Sa,v (constant per cylinder 
Ql< := (i"',^"^^) X K) by the study of functions C7q,d = 4>vPaiPa,v)sa,v piecewise continuous 
in t for all x, constant in x for all volume K, defined as 



N-l 



n=o xer 



One may deduce from the estimates ( |5.5[ ) and (5.13) the following property. 
Lemma 8. (^Space translate of Ua,T>) ■ Under the assumptions (i^|l|) — (-H|6|) . Let V be a 



finite volume discretization of fix (0, T) in the sense of Definition^ and let be a solution 
of ( |3.7| )-( 3rTol ). Then, the following inequality hold: 



(7.1) 



\Ua,v{t,x + y) - Ua,v{t,x)\ dxdt < u:{\y\ 



ln'x{o,T) 

for all y G with il' = {x G Q,, [x, x + y] C fl} and uj{\y\) — >■ when \y\ — t- 0. 
Proof. For a = I and from the definition of Uix), one gets 



/ 

7(0, 



< 



(o.r)xn' 

r 
r 

{0,T)xQ' 

+ 



Ui^vit,x + y) - Ui^v{t,x)\ dxdt 

{pi{pi,v)si,v^ {t, x + y)- {pi{pi,v)si,v^ (t, x) 
sip{t, x + y) (piipLvit, x + y)) - pi{pi^T>it, x)) 



dxdt 
dxdt 



Pi 



i{o,T)xn' 
< El + E2 

where Ei and E2 defined as follows 



{Pi,v){t, x) (si^T>{t, x + y)- si^vit, x) 



dxdt 



(7.2) 



(7.3) 



El = pm 



\si,vit, x + y) - si^T>it, x)\ dxdt, 



Eo 



{o,T)xn' 

Pi{pi,v{t,x + y)) - pi{pi^v{t,x))\ dxdt. 



i{o,T)xn' 

To handle with the space translation on saturation, we use the fact that B^^ is an holder 
function, then 

El < pmC [ \Bisi,vit, x + y))- B{sip{t, x))\^ dxdt 

J{o,T)xn' 

and by application of the Cauchy-Schwarz inequality, we deduce 



Ei<C 



{o,T)xn' 



\B{si^v{t, x + y))- B{sip{t, x))\ dxdt) 
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According to [15j), let y G x G Q' , and L G N{K). We set 



'if IL 



1, if the line segment [x, x + y] intersects a^^i, K and L, 
0, otherwise. 



We observe that (see for more details [15] ) 

(7.4) / liaj,Ax)dx<\aK\L\\y\- 

Jn' 

To simplify the notation, we write ^ instead of 
Now, denote that 

N-l g 

n=0 (7K,L 

N-l g 

<c(\y\Y,^tY, \^k\l\ \b{si,l) - B{si,k)\) ■ 

n=0 crK,L 
1 1 

Let us again write = (d_R',L|o"i^,L|) ^ t^i^^, applying again the Cauchy-Schwarz inequality 



and using the fact that the discrete gradient of the function B is bounded (5.13) to obtain 

(7.5) Ei<C\y\'. 

To treat the space translate of E2, we use the fact that the map p'l is bounded and the 
relationship between the gas pressure and the global pressure, namely : pi = p — p defined in 



(2.5) , then we have 

£;2<max|p^| / \pi^v{t,x + y) -pi^v{t,x)\dxdt 

^ J{o,T)xn' 

(7.6) <max|p^|/ \pTi{t,x + y) - pT){t,x)\dxdt 

^ J{Q,T)xn' 

+ max|p^| / \p{si^v{t,x + y)) - p{si^v{t,x))\dxdt, 

* J{0,T)xQ' 

furthermore one can easily show that p is a C^([0, 1]; M), it follows, there exists a positive 
constant C > such that 



E2<C \pTi{t,x + y) - pD{t,x)\dxdt 

J{o,T)xn' 

+ C \si^x>{t,x + y) - si^T>{t,x)\dxdt. 

J(o,T)xn' 
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The last term in the previous inequahty is proportional to Ei, and consequently it remains 
to show that the space translate on the global pressure is small with y. In fact 

/ \pv{t,x + y)-pv{t,x)\dxdt < V <5t V \pI^'-pI^'\ / /3.,,|Jx)dx 

Jio,T)xn' Jn' 

N-l 

n=0 o"K,L 



Finally, using the fact that the discrete gradient of global pressure is bounded (5.5 ), we deduce 
that 

(7.7) [ \Ui^v{t,x + y)-Ui^v{t,x)\dx<C{\y\ + \yf), 

J{o,T)xn' 

for some constant C > 0. 
In addition, we have 

/ / \Ui,v{t,x + dx)-Ui,v{t,x)\dxdt<2 / \Ui^v{t,x + dx) -Ui,v{t,x)\dxdt 
Jo Jn' Jo Jn' 

+25t [ \U?T.(x)\dx 



where Uj^ = pi{p^)s^ and = {x G dist(x, < \5\}. By (7.7), the assumption 5t — )• as 
size(P) — )• and the boundedness of (JJj^i^)h in L^{i}'g), then the space translates of Ui^v on 
Q' are estimated uniformly for all sequence size(Pm)m tend to zero. 

In the same way, we prove the space translate for a = g. □ 

We state the following lemma on time translate of Ua,v- 

Lemma 9. (^Time translate oftJa,T>)- Under the assumptions {H^ — {H^ . LetV be a finite 
volume discretization of Q x (0,T) in the sense of Definition^ and let pa,T> be a solution of 

(3.7) -(3]l0|. Then, there exists a positive constant C > depending on Q, T such that the 
following inequality hold: 

(7.8) / \tJa,v{t + T,x)-Uc,v{t,x)'\^ dxdt<uj{T), 
Jnx{o,T-T) 

for all T G (0,T). Here Co : — t- M"*" is a modulus of continuity, i.e. limT-^o '^(''") = 0. 

We state without proof the following lemma on time translate of tJap- Following the 
lemma ??, the proof is a direct consequence of the estimations (5.5) and (5.13), then we omit 
it. 

8. Study of the limit 

Proposition 3. Let {T>m)m be a sequence of finite volume discretizations of ft x (0, T) 
such that limm-s.+oo size(Pm.) = 0. Then there exists subsequences, still denoted (so,D„)meN, 
iPa,Vjn)meN verify the following convergences 

(8.1) \\Ua,v^ - Ua,Vm\\L^{n') — > 0, 

(8.2) Ua,Vm — ^ strongly in L^{Qt) and a.e. in Qt for all p > 1, 
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(8.3) Vr,^Bisi,T>J^VB{si) 

(8.4) Vv^pv^^Vp 

(8.5) Sa,V„, > Sa 

(8.6) Pa,Vm >Pa 

Furthermore, 
(8.7) 



B. SAAD AND M. SAAD 

weakly in {L'^{Qt)Y, 
weakly in {L'^{Qt)Y , 
almost everywhere in Qt, 
almost everywhere in Qt- 

< Sq < 1 a.e. in Qt, 

Ua = 4>Pa{Pa)sa a.e. in Qt- 



Proof. For the first convergence (8.1) it is useful to introduce tlie following inequality, for all 

a,b£R, 



[\ea + {l-e)b\d0>l{\a\ + \b\). 
Jo ^ 



Applying this inequality to a = U^^^^ - U^,v^, b 
Uci,Vm, we deduce 

l-T 



u: 



U^T^ , from the definition of 



JQ' 



\Ua,vJt,x) - (t,a;)|dxdt < 2 



T+5t 



lUa^Vmii+^t^x) - Ua,v,^it,x)\dxdt. 



Since 6t tends to zero as size(Pm) — ^ 0, estimate (7.8) in Lemma [o] implies that the right-hand 
side of the above inequality converges to zero as size(Pm) tends to zero, and this established 



.1). 



By the Riesz-Prechet-Kolmogorov compactness criterion, the relative compactness of (C/Q,x)^)meN 
in L^{Qx) is a consequence of the Lemma s [s] and [o} Now, the convergence ( |8.2[ ) in L^{Qx) 
and a.e in Qt becomes a consequence of (8.1). Due to the fact that Ua,Vm is bounded, we 
establish the convergence in L^{Qt). This ensures the following strong convergences 

Pa{Pa,v„,)sa,Vm — ^ L^{Qt) and a.e. in Qt ■ 

Denote by Ua = Pa{Pa)sa- Define the map H : M+ x M+ i-^ M+ x [0,i3(l)] defined by 

(8.9) n{uuug) = {p,B{si)) 

where are solutions of the system 

n{p, B{si)) = pi{p- p{B~\B{si))))B-\B{si)) 
Ug{p,B{si)) = Pa{p-p{B'\B{si))m-B-\B{si)). 

Note that EI is well defined as a diffeomorphism, since 

^ = p',{p-p{B~\B{s,)m-\B{s,))>Q 

= p',{p-p{B-\B{s,m[-p'{B-\B{s,mB-'\B{srmB-\B{s,)) 
+ pi{p-p{B-\B{si))))B-''{B{si))>0 
= -p',ip - p{B-HB{s^)mi - B-\Bis^))) > 



dui 
dB 

dUn 



dp 
dB 



p'2(p-p(^-^(e(^l))))[-;^'(e-l(e(.i)))(fi-l'(fi(.l)))][l-^-l(e(.l))] 

P2{p-p{B~HB{s,mB-'\B{si))<0, 
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and if one of the saturations is zero the other one is one, this conserves that the jacobian 
determinant of the map IHI^"'^ is strictly negative. 



As the map IHI defined in (8.9) is continuous, we deduce 

Pv^ — > P a.e. in Qt, 
S{si^T>m) — > ^* a.e. in Qt- 
Then, as B^^ is continuous, we deduce 

si,Vm — > si = B~^{B*) a.e. in Qt, 



and the convergences (8.5) hold. 

Consequently and due to the relationship between the pressure of each phase and the global 
pressure defined in (2.5), then the convergences (8.6) hold 

Pa.Vm — > Pol a.e. in Qt- 

It follows from Proposition [l] that, the sequence {^VmP'Vm)m&i is bounded in {L'^{Qt)Y, 
and as a consequence of the discrete Poincare inequality, the sequence {pVm)mm is bounded 
in L^{Qt)- Therefore there exist two functions p G LP'{Qt) and V G {L'^{Qt)Y such that 
dsll ) holds and 

^v^Pv^ — > ip weakly in {L'^{Qt)Y. 

It remains to identify Vp by ^p in the sense of distributions. For that, it is enough to show as 
m — )• +oo: 



Em-= / / Vd„^pd^ • 99 dxdt + / / pv^divif dxdt — ^ 0, \/(f e D{Qt). 
J JQt J JQt 

Let T>m be small enough such that 99 vanishes in T^J^ for all K £ T, then 
/ PVmd^y^it,x)dx = ^ / pT>^divip{t,x)dx 

= E Pk [ ^it,x).VK^LdT = ^Yl E iPK-Pl)[ V^(t,x)-r/^|^dr. 

KeT LeN{K) •'''K\L K&TLeN(K) ■'<^k\l 

Now, from the definition of the discrete gradient, 

I ^VmPVmf{t,x)dx =^Y^ I Vv^PVmf{t,x)dx 

= \Y^ E ^(Pl-Pk) [ ip{t,x)-riK\Ldx 

Then, 

-^m = JV V (^k\l{Pl-Pk)( i ^ I / V{t,x) -VKlLdT- ^ I f{t,x) ■'nK\Ldi 

^^TLi]^(K) ^WK\L\Ja^^, \Tk\l\Jt^^^ 
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Due to the smoothness of one gets 

^ ' < C /i, 



f{t,x) ■ ??x|Ldr - I I / ip{t,x) ■ riK\Ldx 



and the Cauchy-Scharwz inequahty with the estimate (5.4) in Proposition [l] yield 

N~l N-1 
\Em\<ChY,^tYl E \cTK\L\\pl-pl<\<ChY,StY, E WK\L\dK\L<Ch\n\T. 

The identification of the hmit in (8.8) follows from the previous convergence. □ 

8.1. Proof of theorem [l| Let T be a fixed positive constant and f G D{[0,T) x fl). Set 
9?^ := for all K £ T and n £ [0,N]. 

For the discrete liquid equation, we multiply the equation (3.8) by 6t(p^^ and sum over 
K €T and n e {0, N}. This yields 

5r + S^ + S^ + ST = 0, 

where 

N-l 

= E E 1^1 {piip^KKK - piip?,KKK) 

n=o xer 

n=0 KeT LeN{K) 
N-l 

n=0 xer LeN{K) 

N-l 

= E E E 1^1 (/'Kpi:;i^)^;:iV]5,ivr' - /'Kp;:i^)(^[,x)'^+v,"ivrO • 

n=0 xer LeN{K) 

Making summation by parts in time and keeping in mind that f(T, xk) = 9^^^^ = 0. For all 
K G T, we get 

N-l 

= - E E 1^1 ^KPi{pi-^')sii' - - E 1^1 '^kpApIk^IkA 

n=o xer i^eTh 
= -EE / / 0i^pKPri^)sri^5t¥'(*'^i^)d2;di - V / 4>kPi{p\^k)^\k'P{^^xk)^x- 

Since (l)v„,Pi{pi.Vm)si,v^ and 4>Vrr^Pi{P^^v^)s^.v^ converge almost everywhere respectively to 
(ppi{pi)si and (t)pi{p^)s^, and as a consequence of Lebesgue dominated convergence theorem, 
we get 

lim Sr= / c^pi{pi)sidMt,x)dxdt- [ (l)pi{p^)s^ip{0,x)dx. 
Now, let us focus on convergence of the degenerate diffusive term to show 
(8.10) lim S^ = - I pi{pi)Mi{si)Vpi-V^dxdt. 

Jq^ 
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Since the discrete gradient of each phase is not bounded, it is not possible to justify the pass 
to the Hmit in a straightforward way. To do this, we use the feature of global pressure and 
the auxiliary pressures defined in (2.5) and the discrete energy estimates in proposition [T] and 
corollary [TJ 

Gathering by edges, the term S"™ can be rewritten as: 

N-l 

n=0 KeT LeN(K) 
N-l 

n=0 KeT LeN{K) 



with, by using the definition (2.5), 
1 

n=0 KeT LeN(K) 



n=0 KeT LeN{K) 



Let us show that 



^.11) lim / piipi)Miisi)Vp-Vipdxdt. 



Qt 

For each couple of neighbours K and L we denote s"^^^ the minimum of and s"^^ and 
we introduce 



n=0 KeT LeN(K) 



Remark that 



N-l 



^r = ^E'^*E E ^\TK\L\PtK]Lms^ 



n+l ..PL -PK f{t'^^^,XL) - f{t'^^^,XK) 



n=0 KeTLeN(K) 
n=0 i^er LeN(K) 

where = 9xk + (1 — 0)xl, < ^ < 1, is some point on the segment ]xk,xl[. Recall that 
the value of ^k\l is directed by r]j^^i, so 

^K\LPV^ ■ VK\L^f{t'''^'^,XK\L) ■ Vk\l = K\LPVm ' V^/J(^+^ 
Define Sa,Vm and s^^j,^ by 
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Now, A^'* can be written under the following continues form 

Jo Jn 

By the monotonicity of B and thanks to the estimate ( 5.13| ), we have 



T 

JQ 



n=0 KeT LeN{K) 
N-1 



ri<\L\ 



n=0 KeT LeN{K) 



d 



K\L 



Bislt') - B{s?+') 



< CsizeiTf. 

Since B^^ is continuous, we deduce up to a subsequence 



.12) 



a.e. on Q^. 



Moreover, we have Sq,-p„ < Sa,Vm ^ ^a,Vm ^^^^ ^a,Vm ~^ a-e. on Qt- Consequently, and 
due to the continuity of the mobility function Mi we have M;(s; Mi{si) a.e on Qt and 

in LP{Qt) for p < +oo. 

As consequence of the convergence (8.6) and by the Lebesgue dominated convergence the- 
orem we get 

Pi{pi,vjMi{s^^j,J{V^)v^ ^ pi{pi)Mi{si)V^ strongly in {L\Qt)Y. 



And as consequence of the weak convergence on global pressure (8.4), we obtain that 



lim A^'*= / pi{pi)Mi{si)Vp-V^dxdt. 

^^+°° Jot 



It remains to show that 
(8.13) 

Remark that 



lim I A 



0. 



< C 



„n+l _ n+1 



Consequently 



< c 



I 

jQ- 



n+l 



Applying the Cauchy-Schwarz inequality, and thanks to the uniform bound on Vv,nPVm and 
the convergence (8.12), we establish (8.13). 



To prove the pass to limit of A^, we need to prove firstly that 

Wminsi)) - ^MK5;;+;j5^2(p(s/))||i2(Q,) ^ O as size(r) ^ 0, 

where r(.0 = /o' yMK^^(z)dz. 

In fact. Remark that there exist a S [si^k^si^l] such as: 
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< C\dl+lipisi))\ < C 



< c 



n+l- 



since B ^ is an Holder function. Thus we get, 



N-1 



n=0 KeT LgN{K) 
N-1 

^E^^E E \TK,L\'-'\TK,Lf\B{slt')-B{sli\ 

n=0 KeTLeN(K) 

and using the Cauchy-Schwarz inequahty and the estimate , we deduce 
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< 



E'^^E E 

n=0 K€T L£N(K) 



1-6 



' AT-l 



n=o _ft:er lgn{k) 



N-l 



=0 KeTLeNiK) ^^^^ 



which shows that ||5^|i'(r(s0) - V^^'(*SL)'^^TL(^(^'))lli2(QT) ^ ^ ^^^^C^) ^ 
from (5.14) in corollary [T| we deduce that there exists a constant C > where the following 
inequalities hold: 

N-l 

E^^E E rK\L{SVil{Tis,))f<C. 

n=0 A'er LeN{K) 



.14) 



V7,„r(s,,2,„J ^ Vr(sO weakly in (L^Qt))^. 



That prove 
(8.15) 
As consequence 

(8.16) ^JMl{sl,vJVv^P{sl,vJ ^ Vr(sz) weakly in (^^(Qr))^ 

Rearranging ^™ to write 

v-i 

= "2 E E E pSil^K^^I JV;,|ip(s;,25„) • VK\LV^{t-+\xK^L) ■ r^K\L, 

n=0 XgT L(^N{K) 

where x^il = ^^j/^- + (1 — 0)xl, < ^ < 1, is some point on the segment ]xk,xl[. using again 
that the mesh is orthogonal, we can write 

Pi{pi.Vm)Mi{sipjVvraP{si,vJ) ' iV(p)v^dxdt. 



Am 
^2 



Jo. 
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As a consequence of the convergences (8.5), (8.6) and by the Lebesgue theorem we get 

Pi iPi,v^ ) \/ Miisi^v^ ) {V^)v^ pi{pi)y/Mi{si)V(p strongly in {L'^{Qt)Y. 
And as consequence of ( |8.16| ), 

(8.17) lim A^ = - r [ piipi)^Mi{si)VT{si) ■ Vifdxdt 
m^+^ Jo Jn 

(8.18) =-/ f pi{pi)Mi{si)Vp{si)-Vipdxdt. 

Jo Jn 

Now, we treat the convergence of the gravity term 

(8.19) lim S^ = - r [ pi{pi)Mi{si)g-Vipdxdt. 

m^ + OO Jq Jq 

Perform integration by parts (3.3) 

N-l 

n=0 KeT LeN{K) 
N-l 

= -2E^*E E F-^J,{^ie-^\x,)-^{e^\xK)). 

n=0 KeT LeN{K) 

Note that the numerical flux F^^^^ is independent of the gradient of pressures and the pass 

to the limit on is mush simple then the term A^'* since the discrete gradient of global 
pressure is replaced by the gravity vector g. We omit this proof of (8.19). 
Finally, can be written equivalently 



^r = EE / / Pl{plt')sl'^'fpit,xMr+\xK)dxdt 

EE/ / Pi{p?i'HKr^'fi{t,xMt-^\xK)dxdt. 



From the convergences (8.5), (8.6) and by the Lebesgue dominated convergence theorem, we 
get 



lim S'^'* = j pi{pi)sifp{t,x)ip{t,x)dx4t pi{pi)sj fi{t,x)(p{t,x)dxdt, 

■m^+oo Jq^ 

which completes the proof of the theorem [l| 



9. Numerical results 

In this section we show some numerical experiments simulating the five spot problem in 
petroleum engineering. A Newton algorithm is implemented to approach the solution of 
nonlinear system (3.8)-(3.9) coupled with a bigradient method to solve linear system arising 
from the Newton algorithm process. 

We will provide two tests made on a nonuniform admissible grid. 
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Figure 2. Mesh with 896 triangles 

Datas used for the numerical tests are the following : 

ki{si) = si, k2{s2) = si 

K = 0.1510-i°m2, (f> = 0.206, 

fi2 = 10^^ Pa.s(water viscosity), fii = 910~^ Pa.s(gas viscosity), 

pip) = Prefi'i- + CrefiP - Pref)), with Pref = 400 Kgm"^, Cref = 10"^Pa~\ Pref = 1.01310^ Pa, 

Lx = Im, Ly = Im (the length and the width of the domain) 

Pc{s) = Pmax{l - S), with Pmax = lO^Pa. 

Initial conditions. Initially the saturation of gas is considered to be equal to 0.9 in the 
whole domain and the gas pressure is considered to be 1.013 10^ Pa. 

Boundary conditions. The wetting fluid (water) is injected in the left-down corner in 
the region ([0, 0.1] x {0}) U ({0} x [0, 0.1]) with a constant pressure equal to 4.026 10^ Pa. The 
right-top corner where ([0.9, 1] x {1}) U ({1} x [0.9, 1]) keeps fluids flow freely at atmospheric 
pressure where as the rest of the boundary is assumed to be impervious (zero fluxes are 
imposed). The influence of boundary conditions can be seen in all figures. 

Meshes. The domain is recovered by 896 admissible triangles see figure [2} 
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Figure 3 . Water field including capillary efi^ect at time T = 6s (left) and at time T = 20s 
with 0.1 < s < 1. 

Figures [3] - 16] show the diffusive effects of the capillary terms, notably the dissipation of 
chocs due to the hyperbolic operator Fig. [6| In fact, during the stage of the displacement 
saturation shock propagate through rock for flows where capillarity terms are neglected, see 
figure [6j This shock, where capillarity effects are signifiant, it is diffused. However, a part of 
the the shock wave maintains its sharp front. 
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